Line width control and trajectory planning for robot guided inkjet deposition

ABSTRACT

In an embodiment, a method for controlling a printer is provided. The method includes: receiving a set of parameters associated with a printer by a computing device, wherein the printer is depositing ink on a substrate to generate a line; measuring a width of the line by the computing device; receiving a duty cycle associated with the printer by the computing device; using the received duty cycle, the set of parameters, and a model, estimating one or more unknown parameters of the set of parameters by the computing device; receiving a desired width of the line by the computing device; and if the desired width is not the same as the measured width: adjusting the duty cycle associated with the printer based on the set of parameters and the model so that the measured width is closer to the desired width by the computing device.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a divisional of U.S. patent Ser. No. 16/913,651, filed on Jun. 26, 2020, entitled “LINE WIDTH CONTROL AND TRAJECTORY PLANNING FOR ROBOT GUIDED INKJET DEPOSITION”, which claims the benefit of priority to U.S. Provisional Patent Application No. 62/868,334, filed on Jun. 28, 2019, entitled “PARAMETER ESTIMATION AND LINE WIDTH CONTROL OF ROBOT GUIDED INKJET DEPOSITION.” The contents of both are hereby incorporated by reference in their entirety.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

This invention was made with government support under award abstract numbers 1933558 and 1563424 awarded by the National Science Foundation. The government has certain rights in the invention.

BACKGROUND

Inkjet deposition (or inkjet printing) refers to multiple technologies to deposit liquids on a substrate. Depositing of various liquid materials has been reported for many applications. Conductive ink has been used to print low-cost antennas on a paper substrate. Organic semiconductors have been printed to create transistors and displays. Inkjet deposition has been used to deposit nerve regeneration conduits, biomaterials and tissues. Solutions of carbon nanotubes have been inkjet printed. A mobile robot equipped with laser scanner and inkjet printer head has been used to autonomously navigate a building and mark structures.

While such applications for inkjet deposition are useful, there is a lack of knowledge regarding how to regulate the placement of a liquid material (e.g., ink) on a substrate (e.g., paper, human tissue, vehicle bodies, etc.) based on changes to drop volume. Current methods for placement regulation rely on different techniques such as position and velocity control, feed-forward control, and meniscus control, each of which is associated with drawbacks.

Furthermore, there are a number of applications that require precise motion planning along a desired trajectory on or over a curved surface, including mobile robots and UAVs conducting coverage path planning, agricultural field automation, vacuum cleaning robots, lawn mowers, etc. There is a particular interest in printing or jetting materials onto curved surfaces, with application in painting robots, printing antennas and printing biological materials. Defining or mapping a specific trajectory to curved surfaces is challenging, and this challenge is increased when the trajectory must pass over specific points and/or remain bounded in a desired region.

SUMMARY

In a first aspect, systems and methods for dynamically controlling the line width of a printer is provided. To achieve this, embodiments may include some or all of the following three components. First, a dynamic model of line width as a function of parameters such as changing duty cycle, nozzle velocity, and the contact angle arising from chemical and physical properties of the liquid (e.g., ink) and surface (e.g., paper) is created. The model may take as an input a current or voltage duty cycle and may output the estimated line width. The model may then be used to obtain control-affine dynamics.

However, as may be appreciated, many of the parameters needed for the model (e.g., terms in the dynamic equation) may be unknown or difficult to measure. As a second component, a nonlinear estimator is used to determine the cumulative uncertainties of the parameters as a single parameter. This single parameter may then be used by the model to estimate line width based on the duty cycle.

The third component is the use of a nonlinear control law to regulate drop diameter. The combination of the three components is a closed-loop regulation of printed line width that operates without any knowledge of the interactions between liquid and surface. This combination of components is robust with respect to parameter error, which guarantees ultimate convergence of error to a known maximum value. In some embodiments, the upper bound for the error may depend on control gain and a difference between real and estimated parameter values. Stability of the estimator and control law may be established through Lyapunov Analysis.

In another aspect, systems and methods for printing a 2D pattern or image on a 3D surface are provided. Part of this approach is the use of constrained optimization to find a surface parametrization that preservers the shape of a desired curve or trajectory from 2D onto a 3D surface. Surface parametrization can be considered as a mapping from a surface to a parametric domain in lower dimensions. Surfaces that are homomorphic to a disk have a parameter space in the 2D plane. This concept lets us map a curve or trajectory defined in 2D to a geometric shape in 3D space. This mapping can result in distortion in the shape and size of the curve. As we want the pattern on the surface to be similar to the original 2D pattern, isometric (length preserving) and conformal (angle preservative) mappings are desirable to minimize distortion.

In an embodiment, a method for controlling a printer is provided. The method includes: receiving a set of parameters associated with a printer by a computing device, wherein the printer is depositing ink on a substrate to generate a line; measuring a width of the line by the computing device; receiving a duty cycle associated with the printer by the computing device; using the received duty cycle, the set of parameters, and a model, estimating one or more unknown parameters of the set of parameters by the computing device; receiving a desired width of the line by the computing device; and if the desired width is not the same as the measured width: adjusting the duty cycle associated with the printer based on the set of parameters and the model so that the measured width is closer to the desired width by the computing device.

In an embodiment, a method for printing a 2D pattern on a 3D surface is provided. The method includes receiving a 2D pattern by a computing device; receiving a point cloud of a 3D surface by the computing device, wherein the point cloud comprises a plurality of points; selecting a subset of the points of the point cloud by the computing device; and generating a mapping of the 2D pattern to the subset of points by the computing device.

This summary is provided to introduce a selection of concepts in a simplified form that are further described below in the detailed description. This summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter.

BRIEF DESCRIPTION OF THE DRAWINGS

The components in the drawings are not necessarily to scale relative to each other. Like reference numerals designate corresponding parts throughout the several views.

FIG. 1 is an illustration of a spot on a level surface;

FIG. 2 is an illustration of a graph;

FIG. 3 is an illustration of a graph;

FIG. 4 is an illustration of a graph;

FIG. 5 is an illustration of a graph;

FIG. 6 is an illustration of a graph;

FIG. 7 is an illustration of an example method for controlling the line width of a printer;

FIG. 8 is an illustration of a discontinuous parameterized curve from 2D to a regular 3D shape;

FIG. 9 is an illustration of a 2D shape mapped to a point cloud corresponding to a sphere;

FIG. 10 is an illustration of a 2D shape mapped to a point cloud corresponding to a manikin cheek;

FIG. 11A is an illustration of a square spiral pattern mapped with angle preservation;

FIG. 11B is an illustration of a square spiral pattern mapped with no angle preservation;

FIG. 12 illustrates the trajectory of the end effector frame tracing a spiral curve on a sphere;

FIG. 13 is an illustration of a method for printing a 2D patterns on a 3D surface; and

FIG. 14 illustrates an example computing device.

DETAILED DESCRIPTION I. Parameter Estimation and Line Width Control of Robot Guided Inkjet Deposition A. Ink-Jet Printing Model Derivation

As described herein, a small amount of liquid in the air having recently left a jetting nozzle of a printer is referred to as a drop, while the liquid in contact with a surface is referred to as a spot. A spot deposited on a solid substrate will flatten and spread. The degree to which a drop of liquid spreads on a solid surface is referred to as wetting. Wetting is a complex phenomenon dependent on the chemical properties of the liquid, solid surface and surrounding air, physical properties of the solid such as surface roughness, temperature, and possible contamination such as dirt, oil, etc. Drop volume has the strongest influence on spot size, followed by fluid viscosity and surface interactions. Drop velocity and impact dynamics are generally negligible for small drops if the impact is normal to the surface. Note while the disclosed systems and methods are described with respect to depositing or printing ink on surfaces, it is for illustrative purposes only. The disclosed systems and methods can be used to print, write, or deposit materials such as liquids or gels on a variety of surfaces.

A second factor influencing line width is drop-to-drop spacing. Constant velocity of the nozzle of the printer mitigates changes in line width due to spacing. Single spot diameter determined by wetting defines the minimum line width, and reduction in spacing between drops/spots increases the line width.

Wetting is typically studied in terms of contact angles (FIG. 1 ), which describe the shape of a spot 101 on a level surface 103. Small contact angles (θ<90°) indicate high wetting and large spread, while large angles indicate low wetting and little spread. Young's equation gives the contact angle θ at equilibrium as a function of the liquid-vapor, solid-vapor, and solid-liquid interfacial tensions (γ_(tv), γ_(sv) and γ_(st) respectively)

${\cos\theta} = {\frac{\gamma_{sv} - \gamma_{sl}}{\gamma_{lv}}.}$

The diameter of a spot on the surface, as viewed from above, can be described as a function of the spot volume V and contact angle. Depending on θ, the diameter d is

$\begin{matrix} {{{{If}\theta} < {90{^\circ}}},{d = {2\left( \frac{3}{\pi} \right)^{1/3}\frac{\sin\theta}{\left( {2 + {\cos\theta}} \right)^{1/3}\left( {1 - {\cos\theta}} \right)^{2/3}}V^{1/3}}}} & (1) \end{matrix}$ ${{{If}\theta} > {90{^\circ}}},{d = {2\left( \frac{3}{\pi} \right)^{1/3}\frac{1}{\left( {2 + {\cos\theta}} \right)^{1/3}\left( {1 - {\cos\theta}} \right)^{2/3}}{V^{1/3}.}}}$

Given accurate knowledge of the material properties, it is possible to estimate θ off-line and determine the correct volume for a drop to give a desired diameter. While research has gone into control to correctly generate the desired volume, if there is error in the estimate of θ, the diameter will be incorrect. Given a desired spot diameter d_(d), the error in spot diameter is defined as e(t)=d(t)−d_(d)(t). Consider the case of θ<90°, and assume θ is constant. The change of spot diameter is a function of the drop volume:

$\begin{matrix} {\overset{.}{d} = {{\frac{dd}{dV}\frac{dV}{dt}} = {\frac{\frac{2}{3}\left( \frac{3}{\pi} \right)^{1/3}\sin\theta V^{- \frac{2}{3}}}{\left( {2 + {\cos\theta}} \right)^{1/3}\left( {1 - {\cos\theta}} \right)^{2/3}}{\overset{.}{V}.}}}} & (2) \end{matrix}$

Collecting the constant terms as

$\begin{matrix} {\beta = \frac{\frac{2}{3}\left( \frac{3}{\pi} \right)^{1/3}\sin\theta}{\left( {2 + {\cos\theta}} \right)^{1/3}\left( {1 - {\cos\theta}} \right)^{2/3}}} & (3) \end{matrix}$

the relationship between {dot over (d)} and {dot over (V)} can be written as

$\begin{matrix} {\overset{.}{d} = {\beta V^{- \frac{2}{3}}{\overset{.}{V}.}}} & (4) \end{matrix}$

Equation (1) also can be represented as

$\begin{matrix} {d = {3\beta V^{\frac{1}{3}}}} & (5) \end{matrix}$ or $\begin{matrix} {V = {\left( \frac{d}{3\beta} \right)^{3}.}} & (6) \end{matrix}$

By substituting (6) in (4) we obtain

$\begin{matrix} {\overset{.}{d} = {9\frac{\beta^{3}}{d^{2}}{\overset{.}{V}.}}} & (7) \end{matrix}$

Consider D(t) as valve duty cycle and estimate the model between {dot over (V)}(t) and {dot over (D)}(t) for a specific range of duty cycle as the linear equation (8). Note that in some embodiments, rather than valve duty cycle voltage, a current duty cycle of a DC motor could be used, as well as proportional current of an AC motor. Equation (8) is

{dot over (V)}(t)=a{dot over (D)}(t)  (8)

in which case {dot over (d)}(t) can be written as

$\begin{matrix} {\overset{.}{d} = {9\frac{\beta^{3}}{d^{2}}\alpha{\overset{.}{D}(t)}}} & (9) \end{matrix}$ or $\begin{matrix} {\overset{.}{d} = {\frac{b}{d^{2}}{\overset{.}{D}(t)}}} & (10) \end{matrix}$

where b=9β³α. To give (10) a control affine form, we design D as

{dot over (D)}=ad ³ +d ² u,a>0  (11)

with u as our new virtual input and a as a parameter that can be chosen by the user to determine the speed of the dynamics. Using (11), we transform (10) to

{dot over (d)}=abd+bu  (12)

Note that b encapsulates all unknown parameters and depends on surface and liquid characteristics. We assume that b is constant for a fixed surface, liquid and valve.

Our final model is then

{dot over (d)}=abd+bu=b(ad+u)=bf(d,u)

γ=d.  (13)

where u can be converted to the real input based on (11).

B. Parameter Estimation and Tracking Control of Drop Spread

In this section, a new, exponentially fast estimation approach is introduced to learn parameters in a class of nonlinear systems. Taking advantage of this knowledge, a Lyapunov-based tracking control is also applied in parallel with the estimation phase to follow the desired value.

1. Preliminaries

Definition 1. A set S⊂

is called measure zero if

∀ϵ>0,∃{I _(k)}_(k=1) ^(∞) ,s.t.S⊆∪ _(k) I _(k)and Σ_(k) l(I _(k))<ϵ,  (14)

where I_(k): =(a_(k), b_(k)) ⊂

are open intervals with length l(I_(k)): =b_(k)−a_(k). In other words, if the sum of the lengths of intervals enclosing all elements of S can be made arbitrarily small, a finite set of points is of measure zero.

2. Parameter Estimation Approach

Theorem 1. For a nonlinear system, we define a dynamic model as

{circumflex over ({dot over (d)})}={circumflex over (b)}f(d,u)  (15)

where {circumflex over (b)}(t)∈

is an estimate of constant b≠0; furthermore, {circumflex over (b)}(t) is updated according to the following rule

{circumflex over ({dot over (b)})}=−k ₁ {circumflex over (b)}+k ₂ s,  (16)

where k₁, k₂∈

₊ are positive scalar gains, and s is given by

$\begin{matrix} {s:=\left\{ \begin{matrix} {\int_{0}^{t}{\left( {\overset{.}{d} - \overset{.}{\hat{d}}} \right){f^{- 1}\left( {d,u} \right)}d\tau}} & {{{if}{f\left( {d,u} \right)}} \neq 0} \\ 0 & {{{if}{f\left( {d,u} \right)}} = 0.} \end{matrix} \right.} & (17) \end{matrix}$

If the set of times

={t∈

|f(d,u)=0} is a measure zero set, then exponentially converges to b from any initial value. Furthermore, if k₁ ²>4k₂, and {circumflex over (b)} has initial condition {circumflex over (b)}(t₀)>0, then {circumflex over (b)} will never equal 0.

Proof. Define our estimation error as

e _(b) =b−{circumflex over (b)}.  (18)

Taking the time derivative of e_(b) gives

ė _(b)={circumflex over ({dot over (b)})}.  (19)

Replacing (16) in (19)

ė _(b) =k ₁ {circumflex over (b)}−k ₂ s.  (20)

Replacing s from (17) when t∉

(i.e., when f (d, u)≠0) gives dynamics

$\begin{matrix} \begin{matrix} {{\overset{.}{e}}_{b} = {{k_{1}\hat{b}} - {k_{2}{\int_{0}^{t}{\left( {\overset{.}{d} - \overset{.}{\hat{d}}} \right){f^{- 1}\left( {d,u} \right)}d\tau}}}}} \\ {= {{k_{1}\hat{b}} - {k_{2}{\int_{0}^{t}{e_{b}{f\left( {d,u} \right)}{f^{- 1}\left( {d,u} \right)}d\tau}}}}} \\ {= {{k_{1}\hat{b}} - {k_{2}{\int_{0}^{t}{e_{b}d{\tau.}}}}}} \end{matrix} & (21) \end{matrix}$

Taking another time derivative of (21) and replacing {circumflex over ({dot over (b)})} with (19) gives

ë _(b) k ₁ ė _(b) k ₂ e _(b)=0.  (22)

Because k₁, k₂>0, then e_(b) exponentially converges to zero and {circumflex over (b)} converges to b.

For t∈

, replacing s from (17) in (16), when f (d, u)=0, the dynamics are given by

{circumflex over ({dot over (b)})}=−k ₁ {circumflex over (b)}  (23)

which guarantees {circumflex over (b)}(t) does not jump or increase ∀t∈

. Considering (18) and that b is constant, e_(b) also does not jump or increase ∀t∈

. As the set of times

is measure zero, those set of times do not contribute to the solution of e_(b)(t), and e_(b) exponentially converges to zero based on (22). Furthermore, if k₁ ²>4k₂, the system described in (22) is overdamped and will not overshoot or experience oscillations. As b>0, {circumflex over (b)} with initial condition {circumflex over (b)}(t₀)>0 will converge to b without crossing 0.

Note that it is possible to estimate d in (17) using a variety of methods. While a backwards difference estimate is generally susceptible to sensor noise, in our case the estimate is immediately integrated, which provides filtering and mitigates noise. Our experimental results use backwards difference and do not show obvious ill effects.

3. Tracking Control Under Stable Parameter Estimation

Let the system in (13) be rearranged to

$\begin{matrix} {\overset{.}{d} = {{{bf}\left( {d,u} \right)} = \begin{matrix} \left\lbrack b \right. & {\left. b \right\rbrack\begin{bmatrix} {ad} \\ u \end{bmatrix}} \end{matrix}}} & (24) \end{matrix}$

Lemma 1. Let d_(d)(t)∈

be a differentiable desired trajectory. Assume b≠0. Under the estimator design given by (16), {circumflex over (b)} is an estimate that exponentially converges to b, according to Theorem 1. Setting u={circumflex over (b)}⁻¹(d_(d)−{circumflex over (b)}ad−ke_(t)), where k∈

⁺ is a positive definite tracking gain. Under the conditions that 1) {circumflex over (b)}(t)≠0, and 2) k₁, k₂ are sufficiently large k₁ ²>>4k₂), then the system in (24) remains stable, and d converges to d_(d) exponentially as time goes to infinity.

Proof. Define e_(t)=d−d_(d), and set a Lyapunov function candidate as

$\begin{matrix} {W = {\frac{1}{2}{e_{t}^{2}.}}} & (25) \end{matrix}$

Taking the time derivative of (25) leads to

$\begin{matrix} \begin{matrix} {\overset{.}{W} = {e_{t}\left( {\overset{.}{d} - {\overset{.}{d}}_{d}} \right)}} \\ {= {e_{t}\left( {{bad} + {b{{\hat{b}}^{- 1}\left( {{\overset{.}{d}}_{d} - {\hat{b}{ad}} - {ke}_{t}} \right)}} - {\overset{.}{d}}_{d}} \right)}} \\ {= {e_{t}\left( {{- b{\hat{b}}^{- 1}{ke}_{t}} + {\left( {{b{\hat{b}}^{- 1}} - 1} \right){\overset{.}{d}}_{d}}} \right)}} \end{matrix} & (26) \end{matrix}$

Defining E_(b)=b/{circumflex over (b)}−1 gives

{dot over (W)}=−k(E _(b)+1)e _(t) ² +E _(b) {dot over (d)} _(d) e _(t)  (27)

By defining C(x, t)=E_(b){dot over (d)}_(d), equation (27) can be rewritten

{dot over (W)}=−k(E _(b)+1)e _(t) ² +C(x,t)e _(t).  (28)

Let i∈

and a_(i), b_(i)∈

⁺. We know from the estimator design and Theorem 1 that e_(b) is exponentially decaying to zero, therefore the absolute value of e_(b) can be upper bounded by a decaying exponential of the form |e_(b)|≤a₀e^(−b) ¹ ^(t). Theorem 1 also showed that {circumflex over (b)}>0 and is bounded from above by some constant. Noting that E_(b)=b/{circumflex over (b)}−1=e_(b)/{circumflex over (b)}, we therefore conclude that |E_(b)|≤a₁e^(−b) ¹ ^(t).

Likewise, C(x, t) can also be upper bounded by an decaying exponential |C(x, t)|≤a₂e^(−b) ² ^(t), and finally equation (28) can be upper bounded as

{dot over (W)}≤(|k|a ₁ e ^(−b) ¹ ^(t) −|k|)|e _(t) ² |+a ₂ e ^(−b) ² ^(t) |e _(t)|

≤−2(|k|−|k|a ₁ e ^(−b) ¹ ^(t))W+√{square root over (2)}a ₂ e ^(−b) ² ^(t)|√{square root over (W)}|  (29)

A closed-form solution for the ODE (29) exists, but it is quite complicated. Instead, we solve W(t) from two different initial conditions, W(t₀)≥1 and W(t₀)<1. We show the solution for the system when W(t₀)≥1 will grow smaller until W(t)<1, and that the solution for the system when W(t₀)<1 grows smaller and is lower bounded by 0. Note then, that if W(t₀)≥1, the system will eventually enter the domain for W(t₀)<1 and thus W(t) converges to 0 for all initial conditions.

For the case W(t)≥1, we have W(t)+√{square root over (W(t))}≤2 W(t), and we can replace (29) with

{dot over (W)}≤−2|k|W+a ₃ e ^(−b) ³ ^(t) W  (30)

where a₃e^(−b) ³ ^(t)>max(2∥K∥a₁e^(−b) ¹ ^(t), √{square root over (2)}a₂e^(−b) ² ^(t)). The ODE (30) with initial condition W(t)=W(t₀) has a solution

$\begin{matrix} {{W(t)} = {{{W\left( t_{0} \right)}\left( e^{- \frac{2a_{3}e^{\frac{2a_{3}}{b_{3}}}}{b_{3}}} \right)e^{- 2{❘k❘}t}} = {{W\left( t_{0} \right)}\gamma e^{- 2{❘k❘}t}}}} & (31) \end{matrix}$

where we replace the constant term in parenthesis with γ. The solution (31) decays exponentially from any W(t₀)>1, and eventually W(t)<1 at some time t₁.

For the case W(t)<1, W(t)+√{square root over (W(t))}≤2√{square root over (W(t))}, and we can replace (29) with

{dot over (W)}≤−2|K|W+a ₃ e ^(−b) ³ ^(t)√{square root over (W)}.  (32)

The ODE (32) with initial condition W(t)=W(t₁) has a solution

$\begin{matrix} {{W(t)} = {\frac{{e^{- 2b_{3}t}\left\lbrack {a_{3} - {\left( {a_{3} - {\left( {b_{3} + {❘k❘}} \right)\sqrt{W\left( t_{1} \right)}}} \right)\left( e^{t({b_{3} - {❘k❘}})} \right)}} \right\rbrack}^{2}}{\left( {b_{3} - {❘k❘}} \right)^{2}}.}} & (33) \end{matrix}$

The solution in (33) decays to zero if |k|≥b₃. Therefore W(t) is upper bounded by an exponentially decaying function for any initial condition W(t₀), and is lower bounded by 0. By the Comparison Lemma [10], W(t)→0 as t→∞, which means d→d_(d).

Note that the requirement that d_(d)(t) be differentiable will rule out some desired patterns, such as piece-wise linear or step changes in line width. These can be approximated as close as necessary with smooth functions, such as splines or smooth heaviside approximations.

Corollary 1

If Lemma 1 is satisfied, and d_(d) does not satisfy

d _(d)(t)=d _(d)(0)e ^(−kt) +e ^(−kt)∫₀ ^(t) kd(s)e ^(−ks) ds.  (34)

for any non-instantaneous interval of time, then the set of times

={t∈

|f(d,u)=0} is a measure zero set for the system described in (13).

Proof. Given f (d, u)=ad+u and replacing u from Lemma 1, we have

f(d,u)=ad+u={circumflex over (b)} ⁻¹({dot over (d)} _(d) k(e _(t))).  (35)

The proof is considered in two cases. In the first case, we assume {dot over (d)}_(d)(t) is not identically equal to zero for a span of time. I.e.,

α≥0, β≤∞,s. t. {dot over (d)}_(d)(t)=0∀t∈[α, β]. Following Lemma 1, e_(t) exponentially converges to zero and {circumflex over (b)}≠0. From Theorem 1, {circumflex over (b)} is bounded and {circumflex over (b)}⁻¹≠0. It remains to show that d _(d)−k(e_(t))=0 only for the set of times

, which is measure 0. Solving the differential equation {dot over (d)}_(d)−k(d−d_(d))=0 leads to (34). As long as d_(d)(t) does not satisfy (34),

t, f (d, u) 0. Note that (34) is a very specific function that is unlikely to occur in practice and can be avoided by the control designer.

In the second case, we assume d_(d)(t) does become identically equal to zero for a span of time. I.e., ∃α≥0, β≤∞,s.t. {dot over (d)}_(d)(t)=0 ∀t∈[α,β]. Because {circumflex over (b)}≠0 and {dot over (d)}_(d)(t)=0, f (d, u) in (35) only depends on e_(t). Note that e_(t) converges to zero exponentially based on Lemma 1. Therefore e_(t) approaches 0 in the limit but is not identically zero for any noninstantaneous interval of time; i.e. for t∈[α,β],

={t|f (d, u)=0} is a measure 0 set.

Theorem 2

Consider a nonlinear system described in (13). We define a dynamic model as (15), where {circumflex over (b)}(t) is the estimate of b≠0 where {circumflex over (b)}(t=0)≥0, and is updated according to (16). k₁ and k₂ are positive and sufficiently large, where k₁ ²>4k₂ designed to render an overdamped second order differential equation of (22), s is defined according to (17), and T is a measure zero set based on Corollary 1. u in (24) can be designed as follows (where we can show {circumflex over (b)}(t) never becomes zero)

u={circumflex over (b)} ⁻¹({dot over (d)} _(d) −{circumflex over (b)}ad−ke _(t)) if {circumflex over (b)}≠0  (36)

Then, {circumflex over (b)} and d remain stable and converge to b and d_(d) respectively, when time goes to infinity.

Proof. Because k₁ ²>4k₂ and b is always positive, {circumflex over (b)}(t) with {circumflex over (b)}(t=0)>0 converges to b in an overdamped fashion and never crosses zero. Substituting (36) and (24) into (26) gives derivative of the Lyapunov function

{dot over (W)}=−2k(E _(b)+1)W+√{square root over (2W)}C(x,t),{circumflex over (b)}≠0.  (37)

This solution obeys the bounds given in (31) for W(t)≥1 and (33) for W(t)≤1. Based on Lemma 1, d remains bounded and converges to d_(d). Based on Theorem 1, and because

in this case is a zero-measure set, {circumflex over (b)} also converges to b independent of the input.

C. Simulation

In this section, we try to simulate and control the line width of the ink jet printing system using our approach and compare our scheme with Model Reference Adaptive Control (MRAC) and Adaptive Pole Placement Control (APPC). Finally, we simulate a set point regulation by applying our approach where MRAC and APPC cannot achieve parameter convergence.

1 Line Width Simulation

In this section, we simulate the proposed estimation and control scheme to regulate line width for both ideal and non-ideal situations.

1.1 Line Width Simulation

We first test our estimator and control scheme under ideal conditions using Simulink. We set b=0.5 and constant. The desired line diameter is a sinusoidal function

d _(d)=sin(t)+1  (38)

which varies between 0 and 2 mm. In the graph 200 of FIG. 2 , one can see that for different initial values of {circumflex over (b)}={0.1, 0.3, 0.7, 0.9}, the estimate converges to the correct value of b. Convergence speed can be tuned through the choice of k₁ and k₂.

In the graph 300 of FIG. 3 , one can see tracking performance of the control law in (36) for b=0.5. We set k=1, a=1. The line 301 is the desired line width, and the line 303 is the real width. The tracking error converges to zero as expected. In the graph 400 of FIG. 4 , the line 401 depicts the control signal. Duty cycle can also be calculated from (11).

1.2 Performance with a Time-Varying Unknown Parameter

In real scenarios, b is not constant due to different factors like dirt, surface non-coherence, etc. We next simulate a situation in which b varies between 0.4 and 0.6 due to the aforementioned factors. At each time instant, we define b as a random variable with mean of 0.5 and variance 0.1 with saturation limits of 0.4 and 0.6. We used the proposed control law with {circumflex over (b)}(t=0)=0.1 and k=1.

The graph 500 of FIG. 5 depicts tracking error with two different tracking gains 501 and 503. For both gains, the tracking errors converge to around zero with little oscillation. As expected, for the higher tracking gain 501, e_(t) converges to zero faster. The graph 600 of FIG. 6 shows the input value 601 (which is approximately sinusoidal, as expected) and line width (desired line width 301 and real line width 303) under parameter uncertainty.

2. Parameter Estimation During Set Point Regulation

Assume the plant in (12). Here, we simulate the plant applying our adaptation law and control method. Set point is r(t)=1. One can observe in contrast to MRAC and APPC, parameter estimation error converges to zero, without PE.

D. Experiments

We present experiments of line width regulation in ink-jet deposition. In one embodiment, a Staubli TX90 robot arm has a solenoid valve is mounted to the end effector. The valve mount was designed in Solidworks and machined such that the offset from the end effector flange is precisely known. We also mounted an Intel RealSense R300 RGB-D camera on the end effector, with image processing routines measuring line width. The robot is controlled to move in a horizontal line with constant desired velocity based on the built-in firmware. Based on the line width error, the proposed estimation in Theorem 1, and control algorithm introduced in Theorem 2, equation (36), are run on C++ code on an Intel i5 laptop. Control outputs are then sent to an Arduino microcontroller, which generates a pulse train to drive the solenoid valve. Colored water under regulated pressure is then deposited. Other computing systems may be used such as the computing system 1400 disclosed in FIG. 14 , for example.

E. Proposed Method

FIG. 7 is an illustration of an example method 700 for controlling the line width of a printer. The method 700 may be implemented by a general purpose computing device in communication with a printer, such as the computing system 1400 described with respect to FIG. 14 . The printer may be an inkjet printer or an extrusion printer. Other types of printers may be supported.

At 710, a set of parameters associated with a printer is received. The parameters may include one or more unknown parameters related to the characteristics of the ink and the surface that is being printed on. Because some of the unknown parameters may be difficult to measure or obtain, a single parameter may be estimated to represent the cumulative uncertainties associated with the set of parameters. The estimated single parameter may be provided by a non-linear estimator.

At 720, a duty cycle associated with the printer is received. The duty cycle may be the valve duty cycle and may represent the amount of time that the valve is open or the speed of a motor. The duty cycle of the printer may be provided by manufacturer of the printer or may be measured by a user or administrator.

At 730, the line width of the line is measured. The line width may be measured by a camera or other sensor capable of measuring the width of a line.

At 740, using the duty cycle, the parameters, and a model, one or more unknown parameters are estimated. Depending on the embodiment, steps 730 and 740 may be performed at or around the same time.

At 750, a desired width of the line is received. The desired with may be provided by a user or administrator. Depending on the embodiment, the desired width may be provided by a particular equation or function that the printer is printing or depositing on the surface.

At 760, it is determined that the measured line width is not the same as the desired line width. The determination may be made by comparing the measured line width to the desired line width. Any method for comparing line widths may be used.

At 770, a duty cycle of the printer is adjusted based on the set of parameters and the model in response to the determination. The duty cycle may be adjusted such that the line width measured by the model is closer to the desired line width. Any method for adjusting the duty cycle of a printer may be used. After adjusting the duty cycle, the method 700 may return to 720 to continue to monitor line width.

II. Trajectory Planning Over Regular General Surfaces with Application in Robot-Guided Deposition Printing

A. Introduction

Systems and methods to design and carry out trajectories over regular curved surfaces are provided. In one embodiment, constrained optimization is used to find a surface parametrization that preservers the shape of a desired curve or trajectory from 2D onto a 3D surface. Surface parametrization can be considered as a mapping from a surface to a parametric domain in lower dimensions. Surfaces that are homomorphic to a disk have a parameter space in the 2D plane. This concept lets us map a curve or trajectory defined in 2D to a geometric shape in 3D space. This mapping can result in distortion in the shape and size of the curve. As we want the pattern on the surface to be similar to the original 2D pattern, isometric (length preserving) and conformal (angle preservative) mappings are desirable to minimize distortion.

The first concept in our approach is to use a Bezier surface parameterization to fit a parametrized surface to a point cloud representation of the surface, such as from a lidar or 3D camera. A Bezier surface is initially defined using points in the point cloud as the control points in the Bezier surface definition. The control points are then adjusted through a novel optimization to provide a better surface fitting between Bezier surface and the point cloud. In this way, there is no need to know the parameterized manifold equation of the surface. The second concept is to extend the surface fitting through a novel constrained optimization procedure such that 2D curves are mapped to the surface with minimal distortions; specifically, the mapping is as close to conformal or isometric as possible. A third concept adds constraints to the optimization that enforce the mapped trajectory lie within a specific region on the surface or pass through specific points. An example of mapping a complicated, discontinuous parameterized curve 810 from 2D to a regular 3D shape 820 is shown in FIG. 8 .

Focusing on our application of deposition printing, we discuss how to take the defined 3D curve and design the trajectory of a robot arm end effector to guide a nozzle along the surface while maintaining a distance to surface and remaining normal to the surface while printing/depositing materials on the curved surface.

B. Mathematical Preliminaries

In this section, we cover a variety of topics from differential geometry and surface fitting that we need in ensuing sections. These topics include regular surfaces, surface parametrization, mapping, conformal mapping, Bezier surfaces, etc.

1. Regular Surfaces and Surface Parametrization

Parametrized surfaces extend the idea of parametrized curves in Euclidean space. It can be considered as a one-to-one mapping from the surface to a specific domain. A parametrized surface is a smooth function r_(s)(u, v): U→

³ where u and v are in

, and U is an open set in

², such that for all q∈U, dr_(s)(q) has rank 2, where dr_(s)(q) is the Jacobian of r_(s) at q.

We are specifically interested in regular surfaces on which the usual notation of calculus on a parametrized surface can be applied. The tangent plane is defined at every point on a regular surface. Consider a surface S⊂

³ with parametric representation

r _(s)(u,v)=[x(u,v),y(u,v),z(u,v)]=r  (1)

for points [u, v] in a topological disk in

² (e.g. a single connected piece with no holes), and r is short notation of the surface parametrization. This representation is regular if for each p∈S, there exist a neighborhood V⊂

³ and a map r_(s):U→V∩S of an open set U⊂

² onto V∩S⊂

³ such that

1. x(u, v), y(u, v), z(u, v) have continuous partial derivatives of all orders in U;

2. r_(s) is a homeomorphism; and

3. For each q E U, the differential dr_(s):

²→

3 is one-to-one.

2. Bezier Surfaces

A regular parametrized surface can be approximated by a two-dimensional Bezier surface. A Bezier surface is a function of two independent variables given by the tensor product equation

P(u,v)=Σ_(j=0) ^(M)Σ_(i=0) ^(N) P _(i,j) B _(i,N)(u)B _(j,m)(v),0≤u,v≤1  (2)

${B_{i,N}(u)} = {\begin{pmatrix} N \\ i \end{pmatrix}{u^{i}\left( {1 - u} \right)}^{N - i}}$

where P_(i,j)∈

³ are control points, and are basis functions (or blending functions) called bi-variate Bernstein polynomials. To derive equation (2) to approximate a given surface, we need (N+1)(M+1) control points chosen from points on the surface. Larger N and M can provide better surface approximation. A Bezier surface will cross through some control points at the corners and lies in the convex hull of all of its control points. Therefore, taking surface points as control points will result in a Bezier surface spline that is not in contact with the surface at the control points. Indeed, the residual error between the Bezier surface spline and 3D surface can be notable. However, we later show that surface points can be an educated initial guess for an optimization problem to fit a Bezier surface to a given surface.

C. Mapping from a Plane onto a Surface

Many pertinent curves/trajectories can be initially designed and defined in 2D and then mapped from 2D onto a given 3D surface through knowledge of a parameterization of the surface. There can be multiple parameterizations for a given surface. We are interested in those that lead to the least possible distortion. For example, when a rectangle is mapped onto a surface, one desires the corner angles remain right. Some distortion indices, such as length, angle and area, are used to represent different mappings.

1. Curve and Surface Parametrization

The image of a map r_(s)(u, v)=[x(u, v), y(u, v), z(u, v)] defines a surface parametrization, where r_(s):

²→

³. The surface is defined by three functions x(u,v), y(u,v),z(u,v); each is a function of two variables u and v. The image of a map r_(c)(θ)=[x(θ), y(θ), z(θ)] from parameter θ∈

to a curve r_(c)⊂

³, is called curve parametrization. This parametrization can be presented in the 2D domain as r_(cp)(θ)=[u(θ), v(θ)], where r_(cp)⊂

² is a curve in a plane such that r_(c)(θ)=r_(s)∘r_(cp)(θ).

As examples of r_(cp) for simple, popular curves, let a, b, ω, α, β ∈

be constants and θ∈[α, β]⊂

.

-   -   Spiral: r_(cp)(θ)=αθ[cos(ωθ), sin(ωθ)]     -   Square spiral:         r_(cp)(θ)=αθ[|cos(ωθ)|cos(ωθ)+|sin(ωθ)|sin(ωθ)|cos(ωθ)|cos(ωθ)−|sin(ωθ)|sin(ωθ)]     -   Boustrophedon path:         r_(cp)(θ)=α[ωt,|cos(ωθ)|cos(ωθ)+|sin(ωθ)|sin(ωθ)]     -   Rose: r_(cp) (θ)=α[cos 2θ cos θ, cos 2θ sin θ].

2. Simultaneous Parametrization and Surface Fitting to Point Clouds

An analytic surface parametrization may be used for mapping a curve onto a surface. Taking advantage of a Bezier surface definition, we present a parametrized function of the surface based on (2), with control points that need to be optimized, such that the optimized Bezier surface represents the given surface. Consider a surface S, for which we seek a parametrization r_(s)(u, v), as in (1). Assume a point cloud has been captured of S, such as from a LiDAR or 3D camera.

-   -   According to the Bezier formulation (with N=M=3 for clarity in         presentation), which can be represented as a tensor product         Bezier patch

$\begin{matrix} {{r_{s}\left( {u,v} \right)} = {{\sum_{j = 0}^{M = 3}{\sum_{i = 0}^{N = 3}{P_{i,j}{B_{i,3}(u)}{B_{j,3}(v)}}}} = \begin{matrix} \left\lbrack 1 \right. & u & u^{2} & {\left. u^{3} \right\rbrack{\underset{K}{\underset{︸}{\begin{bmatrix} K_{0,0} & K_{0,1} & K_{0,2} & K_{0,3} \\ K_{1,0} & K_{1,1} & K_{1,2} & K_{1,3} \\ K_{2,0} & K_{2,1} & K_{2,2} & K_{2,3} \\ K_{3,0} & K_{3,1} & K_{3,2} & K_{3,3} \end{bmatrix}}}\begin{bmatrix} 1 \\ v \\ v^{2} \\ v^{3} \end{bmatrix}}} \end{matrix}}} & (3) \end{matrix}$ where $\begin{matrix} {K = {\begin{bmatrix} 1 & 0 & 0 & 0 \\ {- 3} & 3 & 0 & 0 \\ 3 & {- 6} & 3 & 0 \\ {- 1} & 3 & {- 3} & 1 \end{bmatrix}{\underset{P}{\underset{︸}{\begin{bmatrix} P_{0,0} & P_{0,1} & P_{0,2} & P_{0,3} \\ P_{1,0} & P_{1,1} & P_{1,2} & P_{1,3} \\ P_{2,0} & P_{2,1} & P_{2,2} & P_{2,3} \\ P_{3,0} & P_{3,1} & P_{3,2} & P_{3,3} \end{bmatrix}}}\begin{bmatrix} 1 & {- 3} & 3 & {- 1} \\ 0 & 3 & {- 6} & 3 \\ 0 & 0 & 3 & {- 3} \\ 0 & 0 & 0 & 1 \end{bmatrix}}}} & (4) \end{matrix}$

and P is a tensor of (N+1)(M+1) control points P_(i,j), i∈{0 . . . N}, j∈{0 . . . M}. If P is chosen from points in the point cloud, this Bezier formulation leads to a surface that only passes through some corner points and not all points in the point cloud. This is not a sufficient fit for our purposes. To get a better fit, we find an optimal P, and consequently K, such that the Bezier surface passes through a set of points in the point cloud Pn_(i,j), i∈{0 . . . N}, j∈{0 . . . M}. To realize this optimization, let

${u_{i,j} = \frac{i}{N}},{v_{i,j} = \frac{j}{M}},{i \in \left\{ {0\ldots N} \right\}},{j \in \left\{ {0\ldots M} \right\}}$

be evenly gridded points on the uv plane (note that (N+1)(M+1) is the number of points). Then, let Pn_(i,j) be the point in the point cloud with the projection onto the uv plane closest to the point [u_(i,j), v_(i,j)].

The optimization problem of fitting r_(s)(u_(i,j), v_(i,j)) to Pn_(i,j), i∈{0 . . . N}, j∈{0 . . . M}, can be solved by finding the value of K that minimizes the objective function

$\begin{matrix} {{\underset{K}{\arg\min}J} = {\sum_{j = 0}^{M}{\sum_{i = 0}^{N}{{{{r_{s}\left( {u_{i,j},v_{i,j}} \right)} - {Pn}_{i,j}}}_{2}.}}}} & (5) \end{matrix}$

(recalling that r_(s)(u_(i,j), v_(i,j)) is a function of K in (3) and (4)). To minimize (5), we employ an interior-point optimization method (specifically, we use the fmincon function in MATLAB) with P_(i,j)=Pn_(i,j) as our educated initial estimate. Interior-point optimization is an approach well-known for constrained minimization. This method defines an alternative approximate problem that is easier to solve than the original problem. Although there is no constraint in our optimization problem, we prefer to apply interior-point due to its accuracy and the fact that we will face optimization problems with nonlinear constraints later in this work.

After optimization, we can treat r_(s)(u, v) as a good mapping function from the uv plane onto the surface S, which passes through Pn_(i,j). To achieve still better results, and in case we need to cover a large surface area, we apply piecewise Bezier functions and partition the surface into different subsurfaces

$\begin{matrix} {{r_{s}\left( {u,v} \right)} = \left\{ \begin{matrix} {{r_{s_{1}}\left( {u,v} \right)},{u_{0} \leq u < u_{1}},{v_{0} \leq v < v_{1}}} \\ {{r_{s_{2}}\left( {u,v} \right)},{u_{1} \leq u < u_{2}},{v_{1} \leq v < v_{2}}} \\ \ldots \\ \ldots \\ {{r_{s_{n}}\left( {u,v} \right)},{u_{n - 1} \leq u < u_{n}},{v_{n - 1} \leq v < v_{n}}} \end{matrix} \right.} & (6) \end{matrix}$

and optimization is applied to each subsurface according to (5) plus a continuity condition.

Numerical optimization is sensitive to initial value, and we do not assume convexity of the optimization function. Therefore, the initial guess has a crucial role for fewer number of iterations and to avoid numerical divergence. Using points from the point cloud as the initial control points in Bezier formulation provides an educated initial guess.

3. Conformal Conditions for Fitting Surface

In addition to a good fit, minimum distortion of angles is also desirable. This helps us to keep the structure and geometry of the curve r_(c) similar to the original r_(cp) in 2D. Isometric maps are the ideal, however, in many cases there is no isometric mapping. Conformal maps are a good practical options to be considered. To apply our optimization method to preserve angle, the following statement is presented.

Theorem 1. Let S₁ and S₂ be a pair of regular surfaces. Then, a diffeomorphism r_(s): S₁→S₂ is conformal if it preserves local angles, that is, if for all p∈S₁ and all x, y∈T_(p)S₁, <(x,y)=<(dr_(p)(x),dr_(p)(y)), where dr_(p)(x)=<(∇r)(p),x>, and T_(p)S₁ is the tangent plane of S₁ at the point p.

This leads to a constrained optimization problem

$\begin{matrix} {{\underset{K}{\arg\min}J} = {\sum_{j = 0}^{M}{\sum_{i = 0}^{N}{{{r_{s}\left( {u_{i,j},v_{i,j}} \right)} - {Pn}_{i,j}}}_{2}}}} & (7) \end{matrix}$ ${{s.t}:{\sum_{l = 1}^{L}{{{\angle\left( {x_{l},y_{l}} \right)} - {\angle\left( {{{dr}_{pl}\left( x_{l} \right)},{{dr}_{pl}\left( y_{l} \right)}} \right)}}}_{2}}} = 0$

for some points p_(i) in uv plane (S₁) and x_(l), y_(l)∈T_(p)S₁. This optimization may be solved using an interior-point method, and leads to a mapping fit to points Pn_(i,j) (first term) and preserves local angles based on the second term in (7).

4. Way Points Conditions for Fitting Surface

Having a set of desired way points pt_(i)∈S, i=0, . . . , n, one might want the mapped pattern pass through them. This means, we not only impose that the pattern structure remains similar to the original pattern, but that pattern r_(cp)(θ) should visit the set of points on the surface. This optimization problem can be formulated as

$\begin{matrix} {{\underset{K}{\arg\min}J} = {\sum_{j = 0}^{N}{\sum_{i = 0}^{N}{{{r_{s}\left( {u_{i,j},v_{i,j}} \right)} - {Pn}_{i,j}}}_{2}}}} & (8) \end{matrix}$ ${{s.t}:{\sum_{l = 1}^{L}{{{\angle\left( {x_{l},y_{l}} \right)} - {\angle\left( {{{dr}_{pl}\left( x_{l} \right)},{{dr}_{pl}\left( y_{l} \right)}} \right)}}}_{2}}} = {{0{\sum_{i = 0}^{n}{\min\left( {{{r_{s}\left( {r_{cp}(\theta)} \right)} - {pt}_{i}}}_{2} \right)}}} = 0}$

where K need to be tuned to minimize J for a good fit to the surface points. The second term preserves local angles, and last term forces the mapped pattern to go through the desired points.

D. Experimental Verification of Mapping a Curve to Point Clouds

We present experimental validation of the procedure discussed above.

1. Vision System

In some implementations, a 3D camera such as an Ensenso N35 3D camera by Image Development Systems may be used. The 3D camera uses a stereo-paired of near infrared cameras with an IR projector to generate a dot pattern. The stereo cameras uses the dot pattern to provide texture for matching points in the two images. Both projector and cameras operate at approximately 850 nm wavelength. The result is a 1280×1024 point cloud. Other sized point clouds may be used. At a standoff distance of approximately 0.5 m from a target, the point cloud depths are accurate to about 0.5 mm.

The camera may be calibrated using a calibration board, where the center of the board is considered as the origin of the point clouds generated using the camera after calibration. This coordinate frame is denoted by F_(W), and we assume all point clouds are generated by the camera with respect to (WRT)F_(W).

2 Point Cloud Extraction and Mapping

The above described camera may be used to capture point clouds from a plurality of shapes. In the example described herein the shapes include a sphere and a manikin face. Other shapes may be used.

For each shape, 2500 evenly distributed points Pn_(i,j), i∈{0 . . . 49}, j∈{0 . . . 49} on the surface to cover the desired area were selected. Through these points, an educated initial guess for the control points in Bezier surface is obtained (P_(i), =Pn_(i,j)). Then, by minimizing the objective function in (7), r_(s)(u, v) is achieved as a mapping from 2D plane onto the desired shape. FIG. 9 shows the 2D shape 810 mapped to a point cloud 920 corresponding to a sphere. FIG. 10 shows the 2D shape 810 mapped to a point cloud 1020 corresponding to a manikin cheek. Based on the optimization (7), the face portrait depicted in the 2D shape 310 has a minimal distortion in angles before and after mapping. his is verified, as the angle error metric is less than 0.05 rad in all cases.

To show the necessity of the angle preservation constraint, an example of pattern placement on the manikin cheek is provided. In FIG. 11A, a square spiral pattern 1110 is mapped to a point cloud 1120 with angle preservation based on (7), while FIG. 11B shows the same pattern 1110 mapped to the point cloud 1120 with no angle preservation constraint. As one can see, right angles and the general shape of the pattern 1110 are conserved remarkably better in 11A than in 11B.

E. Path Planning and Guidance of a Robot End Effector for Printing Tasks

In order to use a robot to conduct deposition printing, it may be necessary to use the mapped curve to define a robot trajectory or path. This trajectory must move the print head nozzle along the curve, while maintaining a specific distance to the surface and keeping the nozzle normal to the surface. Similar tasks would be needed for guiding a UAV over a curved surface while keeping sensors aimed at the ground.

In one embodiment, the robot used is a Staubli TX90, with 6 rotational joints. Other robots may be used. A mount for the nozzle was machined, with the offset from the robot flange to the nozzle known through the CAD design. We considered the nozzle tip as the robot end effector. As the camera is calibrated, point clouds are captured WRT board frame (F_(W)). To map trajectories defined in F_(W) to the robot frame F_(R), we use the transformation

$\begin{bmatrix} {\,_{F_{W}}^{F_{R}}R} & {\,_{F_{W}}^{F_{R}}p} \\ 0 & 1 \end{bmatrix},$

where _(F) _(W) ^(F) ^(R) R∈SO(3) and _(F) _(W) ^(F) ^(R) p∈

³ are the orientation and position of F_(W) origin WRT F_(R), respectively.

After coordinate transformation to the robot frame, it is desirable to set the z-axis of the nozzle normal to the surface, and nozzle top position equal to the mapped curve position on the surface plus a z-axis position offset. We define a rotation matrix based on Euler angles as

$\begin{matrix} {{\,_{F_{E}}^{F_{R}}R} = {{{R_{X}(\alpha)}{R_{Y}(\beta)}{R_{Z}(\gamma)}} = \begin{bmatrix} {C{\beta C}\gamma} & {- C{\beta S}\gamma} & {- S\beta} \\ {{C{\alpha S}\gamma} - {C\gamma S\alpha S\beta}} & {{C\alpha C\gamma} + {S\alpha S\beta S\gamma}} & {- S\alpha C\beta} \\ {{S\alpha S\gamma} + {C\alpha S\beta C\gamma}} & {{S\alpha C\gamma} - {C\alpha S\beta S\gamma}} & {C\alpha C\beta} \end{bmatrix}}} & (9) \end{matrix}$

where, C and S stand for cos and sin. Recall that the columns of _(F) _(E) ^(F) ^(R) R are axes of the end effector frame F_(E) in F_(R). The normal vector n∈

³ of the curve on the surface defined in (1) can be calculated as

$\begin{matrix} {a_{n} = \left\lbrack {\frac{\partial{x\left( {u,v} \right)}}{\partial u}\frac{\partial{y\left( {u,v} \right)}}{\partial u}\frac{\partial{z\left( {u,v} \right)}}{\partial u}} \right\rbrack^{T}} & (10) \end{matrix}$ $b_{n} = \left\lbrack {\frac{\partial{x\left( {u,v} \right)}}{\partial v}\frac{\partial{y\left( {u,v} \right)}}{\partial v}\frac{\partial{z\left( {u,v} \right)}}{\partial v}} \right\rbrack^{T}$ $n = \frac{a_{n} \times b_{n}}{\sqrt{{❘a_{n}❘}^{2}{❘b_{n}❘}^{2}} - {❘{a_{n} \cdot b_{n}}❘}^{2}}$

where u and v are parameterized as [u, v]=[u(θ), v(θ)] for θ⊂

. By setting the last column of the rotation matrix in (9) to n, we set the z direction of the end effector equal to surface normal vector; α and β are then derived as

α=atan 2(−n ₂ ,n ₃),β=atan 2(n ₁ sin(α),n ₂)  (11)

As the end effector moves along the curve, we need to avoid extreme rotation of the final robot joint to prevent tangling/pulling the tubes leading to the nozzle and prevent the robot from hitting its joint limit. We achieve this by setting the second element of the first column of the rotation matrix in (9) (i.e., the x-axis of the end effector WRT F_(R)) equal to zero. This enforces the end effector x-axis to remain tangent to the xz plane in F_(R). This leads to

γ=atan 2(sin(α)sin(β), cos(γ)).  (12)

Given α, β, γ and n, the axes F_(E) are uniquely defined for every point on the curve.

FIG. 12 illustrates the trajectory of the end effector frame tracing a spiral curve on a sphere, where the axes 1201 is the z-axis along the surface normal and the axis 1203 is the x-axis, which always stays aligned with the same plane. Knowing the rotation (R) and position (r_(s)(u, v)=p with an offset along the surface normal vector), one knows desired robot end effector position and orientation to draw the desired path.

F. Proposed Method

FIG. 13 is an illustration of a method 1300 for printing a 2D patterns on a 3D surface. The method 1300 may be performed by one or more general purpose computing devices such as the computing system 1400.

At 1310, a 2D pattern is received. The 2D pattern may be an image and may be received from a user or administrator. Any type of 2D pattern may be supported.

At 1320, a point cloud of a 3D surface is received. The point cloud may be received from a 3D camera and may be a 1024×1024 point cloud. Other types of point clouds may be supported. The point cloud may have a coordinate frame that is denoted as F_(W).

At 1330, a subset of the points of the point cloud are selected. The subset of points are the initial control points and may be selected as described above for placing the 2D pattern on the 3D surface.

At 1340, a mapping of the 2D pattern to the subset of points is generated. The mapping r_(s)(u, v) of the 2D pattern may be mapping from the 2D pattern to the points of the control points corresponding to the 3D surface. The mapping may be generated by minimizing the objective function described at (7). Other methods may be used. The mapping may be both isometric and conformal.

At 1350, a path is generated on the 3D surface using the mapping. The path may be a trajectory that will be followed by the print head nozzle of the printer over the 3D surface while maintaining a specified distance from the 3D surface. The path or trajectory may be determined as described above. In some embodiments, the path may be defined as a .csv file that lists the coordinate of the Cartesian frame [(x,y,z,α,β,γ)] at each time step. Other types of file or definitions may be used.

At 1360, a printer is instructed to draw the 2D pattern on the 3D surface using the generated path. Any method for providing instructions to a printer or robot may be used.

FIG. 14 shows an exemplary computing environment in which example embodiments and aspects may be implemented. The computing device environment is only one example of a suitable computing environment and is not intended to suggest any limitation as to the scope of use or functionality.

Numerous other general purpose or special purpose computing devices environments or configurations may be used. Examples of well-known computing devices, environments, and/or configurations that may be suitable for use include, but are not limited to, personal computers, server computers, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, network personal computers (PCs), minicomputers, mainframe computers, embedded systems, distributed computing environments that include any of the above systems or devices, and the like.

Computer-executable instructions, such as program modules, being executed by a computer may be used. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. Distributed computing environments may be used where tasks are performed by remote processing devices that are linked through a communications network or other data transmission medium. In a distributed computing environment, program modules and other data may be located in both local and remote computer storage media including memory storage devices.

With reference to FIG. 14 , an exemplary system for implementing aspects described herein includes a computing device, such as computing device 1400. In its most basic configuration, computing device 1400 typically includes at least one processing unit 1402 and memory 1404. Depending on the exact configuration and type of computing device, memory 1404 may be volatile (such as random access memory (RAM)), non-volatile (such as read-only memory (ROM), flash memory, etc.), or some combination of the two. This most basic configuration is illustrated in FIG. 14 by dashed line 1406.

Computing device 1400 may have additional features/functionality. For example, computing device 1400 may include additional storage (removable and/or non-removable) including, but not limited to, magnetic or optical disks or tape. Such additional storage is illustrated in FIG. 14 by removable storage 1408 and non-removable storage 1410.

Computing device 1400 typically includes a variety of computer readable media. Computer readable media can be any available media that can be accessed by the device 600 and includes both volatile and non-volatile media, removable and non-removable media.

Computer storage media include volatile and non-volatile, and removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Memory 1404, removable storage 1408, and non-removable storage 1410 are all examples of computer storage media. Computer storage media include, but are not limited to, RAM, ROM, electrically erasable program read-only memory (EEPROM), flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by computing device 900. Any such computer storage media may be part of computing device 1400.

Computing device 1400 may contain communication connection(s) 1412 that allow the device to communicate with other devices. Computing device 1400 may also have input device(s) 1414 such as a keyboard, mouse, pen, voice input device, touch input device, etc. Output device(s) 1416 such as a display, speakers, printer, etc. may also be included. All these devices are well known in the art and need not be discussed at length here.

It should be understood that the various techniques described herein may be implemented in connection with hardware components or software components or, where appropriate, with a combination of both. Illustrative types of hardware components that can be used include Field-programmable Gate Arrays (FPGAs), Application-specific Integrated Circuits (ASICs), Application-specific Standard Products (ASSPs), System-on-a-chip systems (SOCs), Complex Programmable Logic Devices (CPLDs), etc. The methods and apparatus of the presently disclosed subject matter, or certain aspects or portions thereof, may take the form of program code (i.e., instructions) embodied in tangible media, such as floppy diskettes, CD-ROMs, hard drives, or any other machine-readable storage medium where, when the program code is loaded into and executed by a machine, such as a computer, the machine becomes an apparatus for practicing the presently disclosed subject matter.

Although exemplary implementations may refer to utilizing aspects of the presently disclosed subject matter in the context of one or more stand-alone computer systems, the subject matter is not so limited, but rather may be implemented in connection with any computing environment, such as a network or distributed computing environment. Still further, aspects of the presently disclosed subject matter may be implemented in or across a plurality of processing chips or devices, and storage may similarly be effected across a plurality of devices. Such devices might include personal computers, network servers, and handheld devices, for example.

Although the subject matter has been described in language specific to structural features and/or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the specific features or acts described above. Rather, the specific features and acts described above are disclosed as example forms of implementing the claims. 

What is claimed:
 1. A method for printing a 2D pattern on a 3D surface comprising: receiving a 2D pattern by a computing device; receiving a point cloud of a 3D surface by the computing device, wherein the point cloud comprises a plurality of points; selecting a subset of the points of the point cloud by the computing device; and generating a mapping of the 2D pattern to the subset of points by the computing device.
 2. The method of claim 1, wherein generating the mapping of the 2D pattern to subset of points comprises solving an optimization problem.
 3. The method of claim 1, wherein the mapping is isometric and conformal.
 4. The method of claim 1, wherein generating the mapping of the 2D pattern to subset of points comprises solving the optimization problem: ${\underset{K}{\arg\min}J} = {\sum\limits_{j = 0}^{M}{\sum\limits_{i = 0}^{N}{{{r_{s}\left( {u_{i,j},v_{i,j}} \right)} - {Pn}_{i,j}}}_{2}}}$ ${{s.t}:{\sum\limits_{l = 1}^{L}{{{\angle\left( {x_{l},y_{l}} \right)} - {\angle\left( {{{dr}_{pl}\left( x_{l} \right)},{{dr}_{pl}\left( y_{l} \right)}} \right)}}}_{2}}} = 0$
 5. The method of claim 1, further comprising: generating a path on the 3D surface using generated mapping; and instructing a printer to draw the 2D pattern on the 3D surface using the generated path.
 6. The method of claim 5, wherein the path comprises a trajectory that to be followed by a print head nozzle of a printer over the 3D surface while the print head nozzle maintains a specified distance from the 3D surface.
 7. The method of claim 1, wherein the path comprises coordinate of a Cartesian frame [(x,y,z,α,β,γ)] at each time step of a plurality of time steps.
 8. A system for printing a 2D pattern on a 3D surface comprising: a printer; and a computing device adapted to: receive a 2D pattern; receive a point cloud of a 3D surface, wherein the point cloud comprises a plurality of points; select a subset of the points of the point cloud; and generate a mapping of the 2D pattern to the subset of points.
 9. The system of claim 8, wherein generating the mapping of the 2D pattern to subset of points comprises solving an optimization problem.
 10. The system of claim 8, wherein the mapping is isometric and conformal.
 11. The system of claim 8, wherein generating the mapping of the 2D pattern to subset of points comprises solving the optimization problem: ${\underset{K}{\arg\min}J} = {\sum\limits_{j = 0}^{M}{\sum\limits_{i = 0}^{N}{{{r_{s}\left( {u_{i,j},v_{i,j}} \right)} - {Pn}_{i,j}}}_{2}}}$ ${{s.t}:{\sum\limits_{l = 1}^{L}{{{\angle\left( {x_{l},y_{l}} \right)} - {\angle\left( {{{dr}_{pl}\left( x_{l} \right)},{{dr}_{pl}\left( y_{l} \right)}} \right)}}}_{2}}} = 0$
 12. The system of claim 8, wherein the computing device is further adapted to: generate a path on the 3D surface using generated mapping; and instruct a printer to draw the 2D pattern on the 3D surface using the generated path.
 13. The system of claim 12, wherein the path comprises a trajectory that to be followed by a print head nozzle of the printer over the 3D surface while the print head nozzle maintains a specified distance from the 3D surface.
 14. The system of claim 8, wherein the path comprises coordinate of a Cartesian frame [(x,y,z,α,β,γ)] at each time step of a plurality of time steps.
 15. A computer-readable medium with computer-executable instructions stored thereon that when executed by a computing device cause the computing device to: receive a 2D pattern; receive a point cloud of a 3D surface, wherein the point cloud comprises a plurality of points; select a subset of the points of the point cloud; and generate a mapping of the 2D pattern to the subset of points.
 16. The computer-readable medium of claim 15, wherein generating the mapping of the 2D pattern to subset of points comprises solving an optimization problem.
 17. The computer-readable medium of claim 15, wherein the mapping is isometric and conformal.
 18. The computer-readable medium of claim 15, wherein generating the mapping of the 2D pattern to subset of points comprises solving the optimization problem: ${\underset{K}{\arg\min}J} = {\sum\limits_{j = 0}^{M}{\sum\limits_{i = 0}^{N}{{{r_{s}\left( {u_{i,j},v_{i,j}} \right)} - {Pn}_{i,j}}}_{2}}}$ ${{s.t}:{\sum\limits_{l = 1}^{L}{{{\angle\left( {x_{l},y_{l}} \right)} - {\angle\left( {{{dr}_{pl}\left( x_{l} \right)},{{dr}_{pl}\left( y_{l} \right)}} \right)}}}_{2}}} = 0$
 19. The computer-readable medium of claim 15, wherein the computing device is further adapted to: generate a path on the 3D surface using generated mapping; and instruct a printer to draw the 2D pattern on the 3D surface using the generated path.
 20. The computer-readable medium of claim 19, wherein the path comprises a trajectory that to be followed by a print head nozzle of the printer over the 3D surface while the print head nozzle maintains a specified distance from the 3D surface. 